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Abstract 

Thin, penetrable electromagnetic inclusions buried within a half space are retrieved from 
Multi-Static Response (MSR) matrix data collected above this half space at several fre- 
quencies. A non-iterative algorithm is proposed to that aim. It relies on the modeling of 
the MSR matrix according to a rigorous asymptotic expansion of the scattering amplitude. 
The supporting curves of single as well as multiple thin inclusions are retrieved in fairly 
good fashion. This is illustrated by a set of numerical inversions that are carried out from 
synthetic data, both noiseless and noisy ones. 
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1. Introduction 

To achieve a reliable imaging of thin inclusions and cracks still remains a challenging 
problem in the area of electromagnetic non-destructive evaluation. This is true even if 
simple models of those defects are chosen, e.g., each is assimilated with a thin (with respect 
to the wave length of the probing electromagnetic field), homogeneous penetrable layer which 
extends around a properly smooth supporting curve, transmission conditions being satisfied 
at its boundary, with the limiting case of impenetrability being taken care of as well. Yet, 
ill-posedness and inherent nonlinearity of the inverse problem at hand often impair that 
achievement. 

So, over a number of decades, many imaging methods have been suggested. Almost 
all are based on least-square minimization iterative schemes, so consequently, initial guess 
close enough to the unknown defect and suitable regularization tuned to the problem at 
hand are essential. Without, one might suffer from the occurrence of several minima, large 
computational costs, and so on. 



* Corresponding author. 
Email addresses: parkwk@kookmin.ac.kr (Won-Kwang Park), 
dominique . lesselierOlss . supelec . f r (Dominique Lesselier) 

Preprint submitted to Elsevier November 6, 2012 



Recently, MUSIC (Multiple Signal Classification)-type, non-iterative algorithms have 
been developed in order to overcome some of these difficulties. They aim at the retrieval of 
the location and shape of penetrable, electromagnetically thin inclusions as well as the one 



of perfectly conducting ones at a fixed single frequency. From several results in [23l . \24l |25 
such MUSIC algorithms appear both fast and robust, while easily extending to several 
disjoint inclusions. However, whenever applied to limited view data (one cannot illuminate 
the unknown inclusions from all around neither collect the scattered fields all around), for 
example, when those inclusions are embedded within a ha lf sp ace, emitters and receivers 



being set above it, retrievals are much poorer, as is seen in jsl, [21 , 22 , 25 ] for either volumetric 
inclusions or thin ones. Yet, recently, an effective and robust multi- frequency algorithm has 
been proposed so as to find the location of small cracks j3] but a suitable algorithm for 
extended ones appears still in need. 

The purpose of the present paper is to propose an effective, non-iterative imaging algo- 
rithm which is able to work on limited view data in order to retrieve electromagnetically 
penetrable thin inclusions buried within a homogeneous (lower) half space. It is based on 
the fact that the Multi-Static Response (MSR) matrix which one can collect in that set- 
ting (at least samples of it) can be modeled via a rigorously derived asymptotic expansion 
formula of the scattering amplitude in the presence of the inclusions. A number of numeri- 
cal simulations will then illustrate how the proposed imaging algorithm operated at several 
frequencies behaves, and enhances theperformance of imaging that would be based on the 
statistical hypothesis testing, refer to |4|. 

The paper is structured as follows. In section |2j the two-dimensional direct scattering 
problem is sketched and the asymptotic formula for the scattering amplitude are introduced. 
In section [3j the multi-frequency based non-iterative imaging algorithm is outlined. In 
section HI the numerical results are shown. A short conclusion follows (section [5]) . 

2. Direct scattering problem 

In this section, we briefly discuss the two-dimensional, time-harmonic electromagnetic 
scattering from a thin inclusion buried in a homogeneous half space. A more detailed 
description is found in 

Let us define the lower half space and the upper half space represented as 

R 2 _ = {( Xl ,x 2 ) T e M 2 : x 2 < 0} and R 2 + = {(x u x 2 ) T G M 2 : x 2 > 0} , 
respectively. 

A thin penetrable inclusion, T is fully embedded in the homogeneous space M 2 ,. It is 
curve-like, i.e., it is localized in the neighborhood of a curve: 

r = {x + r]n(x) : x G a, 77 G (—h, h)} , 

where a is a simple, smooth curve contained in R^_, with strictly positive distance from the 

boundary, dM. 2 = {(xi,x 2 ) T 6 I 2 : i 2 = 0}, n(x) is a unit normal to a at x, and h is a 

small (with respect to the wavelength of the electromagnetic field in the embedding space 
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Figure 1: Sketch of the thin inclusion T buried within a half-space 



at the given frequency of operation u), positive constant which specifies the thickness of the 
inclusion, refer to Fig. [TJ 

All materials are characterized by their dielectric permittivity and magnetic permeability 
at uj; e_, e + and Et denote the electric permittivity of M 2 ., M.+ and T, respectively. The mag- 
netic permeability can be denoted likewise. So, the piecewise-constant electric permittivity 
< e(x) < +00 and magnetic permeability < < +00 as 

{£+ for x G IR+ ( /i + for x G R 2 + 

e_ for x G 1R 2 \r and fj,(x) = I [i- for x G IR 2 \f 
St for x G Y y \i T for x G T. 

For convenience, the electric permittivity Sq(x) and magnetic permeability fio(x) when there 
are no inclusions are such as 

/ \ f e + for x G Ei , , s f u+ for x G IRl 

£o(:E) = ( e_ for x G ^ = U- for x G M 2 ^ 

Accordingly, the piecewise positive real- valued wavenumber k(x) when there are no inclu- 
sions can be written as 



k(x) 



k + = Uy/e + fi + for x G 
k_ = u^/eZJiZ for x G M 2 ,. 



Let now # = (Q\, 62) be a two-dimensional vector on the unit circle S l C M 2 and Ui(x) = 
e ik + e-x ^ e a p] anar incident wavefield generated in the upper half space on At u, let 
denote the total field which satisfies the Helmholtz equation 

V ■ ( . . u(x)\ + u 2 e(x)u(x) =0 in 



ji(x) 

transmission conditions holding at boundaries cM 2 and dT. 

Let Uq(x) be the plane- wave solution to the Helmholtz equation in the absence of in- 
clusions. Then, within each half space, u s (x) = u(x) — uo(x) is required to satisfy the 
Sommerfeld radiation condition 



i- a — r ( du s (x) . . . A 

hm v \x\ — —, — ; ik{x)u s {x) = 

|*|->oo v 1 1 V d\x\ v ; v 7 



uniformly into all directions x = #r. 

The scattering amplitude is defined as a function K(y, 6) that satisfies 

e *%l 

u s {y) = ^=K(y,6) + o 



y\ \v\y 

as \y\ — > oo uniformly on y = -A and 9 G S 1 . 

In order to express the asymptotic formula for K(y, 6), let us denote v G C 2 and T(x) G C 

as 

2/z_££ 2 

/z_£x 2 + /i+sign(x 2 )A/ 1 - £ 2 ^i 



v(z) = U£i,sign(x 2 Wl - £ 2 x 2 , T(x) 



and a positive, symmetric matrix A(x) satisfying (see [8|, [HJ, 12. 

• *4.(x) has eigenvectors t(x) and n(x) 

• The eigenvalue corresponding to t(x) is 2/^_ 



• The eigenvalue corresponding to n(x) is 2/x_ 

where t(x) and are unit vectors respectively tangential and normal to x G o~. 

Let us emphasize that in the problem at hand, the measurement of incident and observed 
fields is restricted to the upper half space R 2 + . For convenience, we divide the unit circle S 1 
into 

S\ = {x G R 2 + : \x\ = 1} and Si = {x G R 2 _ : \x\ = 1} . 

In that configuration, and using material from 0, 0, [H], EH 12 , we obtain the following 
result. 

Theorem 2.1. For every y G S]_ and 9 £ Si, the asymptotic formula for the scattering 
amplitude K(y, 8) is expressed as 



K%6) =h ± - J ny)T(9) 



v{y) ■ A{x) ■ v(9)e- ik - (v ^- vm - x da(x) 
where the remaining term o(h) is independent of points x G a. 



(1) 

+ o(h), 



Let us notice that the asymptotic expansion (jTJ) is not obtained by mere application of 
the Born approximation. One should refer to [6| for a similar analysis in the case of small 
volumetric inclusions. 
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3. Non-iterative imaging algorithm 



In this section, we apply the asymptotic formula for the scattering amplitude ([T]) so as 
to build up the imaging algorithm. To do so, we use the eigenvalue structure of the Multi- 
Static Response (MSR) matrix K = (Kji), whose element Kji is the amplitude collected at 
observation number j for the incident wave numbered /. 

Let us assume that for a given frequency u, the thin inclusion is divided into M different 
segments of size of order ~. Having in mind the Rayleigh resolution limit from far- field data, 
any detail less than one-half of the wavelength cannot be retrieved, and only one point, say 
x m for m — 1, 2, - • • , M, at each segment is expected to contribute to the image space of 
the response matrix /C, refer to 0, US 23, 24, 25 



For simplicity, let us denote the constant C = h 



(fc_) 2 M+ (i+0 

4/i_ \l fc_|_7T 



and remove the residue term 



o(h) from ([I]). Then, for each y~j = —9j, the jl— th element of the MSR matrix Kji G C, 
j,l = 1,2,--- ,N, is 



K 5l =CK{y v 9 l 
=CT(^)T(^) 



o ik-(v(ej)+v(0i))-x 



da(x) 



M 



m=l L 



fl -lj +2^ -l)/^,)-r(.r () ,)/^,) •-(,•„ 



[Mr 



3 ifc_(«(fl 3 -)+v(tf,)) 



where |cr| denotes the length of a. 

Notice that MSR matrix /C can be decomposed as follows: 

K, = VSV T . 

Here, the matrix S G M. 3Mx3M is a diagonal matrix with component 

(/c_) 2 /i+(l + i) W\ ( E £ 



£ = h- 



where 



E e = M x M diagonal matrix with components 



Ea = 2M x 2M diagonal matrix with 2x2 blocks 



and the matrix T> G C 7Vx3M can be written as follows 



- 1 



2(^-1 



2(1-^ 



[Dl Dl -.. D? Dl Dl 
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Dl M ] 



(2) 



(3) 



where 

D™ = (T(6 1 )e ik - vi6l > Xm ,T(e 2 )e ik - v( - e2) - Xm , ■■■ , T(6 N )e ik - v{eNyXm ) T , 
D 2(m-i)+. = ^ s{Xm) . v^T^e*-^*™, • • ■ , Ux m ) ■ v(9 N )T(6 N )e lk - v ^> x ™) T 

with 

J r(x m ) if s = 1 
U^J .- | nW , f s = 2 

Let us notice that with the representation ([3]), /C is symmetric but not Hermitian (a 
Hermitian matrix could be formed as /C/C). Since /C is is not self-adjoint, a Singular Value 
Decomposition (SVD) has to be used instead of the eigenvalue decomposition. Let us per- 
form this decomposition of matrix /C and let M be the number of nonzero singular values 
for the given u. Then, K can be represented as follows: 

M 

K, =U{u)S{u)V{u) « ^ «m(w)Sm(w)C(w), 

m=l 

where s m (u;) are the singular values, u m (u) and v m (u) are the left and right singular vectors 
of K for m = 1,2,--- ,M. 

Based on the above singular value decomposition, the imaging algorithm is developed as 
follows. For c G M 3 \ {0}, let us define a vector 

u) = (c • (1, uCflOWOe*-**^ ■ • • , c • (1, ^^mM^^f 
and corresponding normalized vector o>) = pjf^yu • Then for m = 1, 2, • • • , M, 



u m ~ e 111 d(x m ; u) and v m ~ e tl2 d(x m ; u) (4) 



with 71 + 72 = f, refer to |l8j ]. Since the first M columns of the matrix U(u) and V(cu), 
{ui(u)), u 2 (oo), ■ ■ • ,um(u)} and {vi(oo), v 2 (co), ■ • • ,%(tt))}, are orthonormal, one can easily 
find that 

(d(x;u),u m (u))) ^ and (d(x;u)),v m (u)) / if x = x m 

(5) 

w), it m (o;)) ~ and (d(x; u), v m (u)) ~ if i/i m 

for m = 1, 2, ■ • • , M, where (a, b) = a ■ b. 

Now, let us construct the following image function with MSR matrix at fixed frequency 

A/ 

W ^( a; ) = w )^m(w))(rf(x;a;),U m (a;))|. (6) 

771=1 

Based on 03]) and (J5J), the map of Wg(x) is expected to exhibit peaks of magnitude of 1 at 

location x m for m = 1, 2, • • • , M and of small magnitude at x G M?_\r. 
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Unfortunately, the image functional fl6]) at single frequency offers an image with poor 



resolution, refer to |4|, [21|, |22j. In order to improve the imaging performance, we suggest a 



normalized image functional at several frequencies {ojj : / = 1, 2, • • • , F}: 

W F (x) = — \(d(x;u f ),u m (u f )}(d(x;u f ),v m (u f )}\, (7) 
f=l m=l 

where Mf is the number of nonzero singular values of the MSR matrix at for / = 
1,2, - ■ ■ , F. Then, similarly with the flS]), the map of W(x) should exhibit peaks of mag- 
nitude of 1 at location x m for m = 1,2,- •• ,M/, and of small magnitude at x G M?_\r. 
A suitable number of Mf for each frequency Uf can be found via careful thresholding, see 



23l l24j . |25| for instance. It is worth mentioning that multiple frequencies should enhance 
the imaging performance via higher signal-to-noise ratio (SNR) based on the statistical 
hypothesis testing, refer to appendix ?? (see Q also). 

Remark 3.1. For effective imaging, to properly choose the vector c = (0,61,62) G M 3 \ {0} 
is a strong prerequisite. The choice comes from the structure of the elements of the MSR 
matrix (TJ|). First, c = (1,0,0) is a good, easy choice for a purely dielectric contrast case. 
But, for a purely magnetic contrast case, the vector b must be a linear combination of the 
tangential and normal vectors while keeping a = 0. To appraise the tangential (and also 
normal) vector of the supporting curve a should not be straightforward, yet we can obtain a 
good image with (61, 62) = (1,0) and (61, 62) = (0, 1) for /it < and /iy > respectively. 
For a double contrast case, c should be a combination of the above two results. A detailed 
discussion can be found in plqL section 4-3.1]. 



4. Numerical examples 

In this section, numerical examples are provided. Throughout, the thickness h of the thin 
inclusion are set to 0.015 and the applied frequency is Uf — jj] here A/, / = 1, 2, • • • , F, is 
the given wavelength. In this paper, frequencies Uf are equi-distributed within the interval 
As for the observation directions yj, they are taken as 

Vj = - (cos Cj, sin Q) , Q = a + (f3 - oc)j^—. -, 

where a = | and (3 = for j = 1,2, • • • , N. See Fig. [2] for an illustration of the test 
configuration. 

Note that when the upper half space is more refractive than the lower one, i.e., k + > k_, 
the number N + of propagating transmitted waves might be less than N, refer to Table [1] 
(see jof also). 

Two CTj characteristic of the thin inclusion Tj are chosen for illustration: 

o~ 1 = {(z-0.2,-0.52 2 -1.5) : zE [-0.5,0.5]) 
a 2 = {(z + 0.2,2 3 + 2 2 -2.5) : z G [-0.5,0.5]} 

7 



Figure 2: Sketch of the test configuration. 



inclusion 


N 


N + 


F 


frequency range 


search domain 


r 2 

r 

L M 


32 
40 
48 


24 
28 
32 


30 
36 
40 


w G [Sj p] 


n = [-i,i] x [-3,-1] 
n = [-i,i] x [-3,-1] 
n = [-1,1] x [-3,-1] 



Table f : (Permittivity contrast case) Test configuration for Ti, T2 and T M . 



and define r M = Ti U T 2 for multiple inclusions. Throughout this section, we denote Sj and 
be the permittivities and permeabilities of Tj, respectively. 

It is worth emphasizing that the data set of the MSR matrix K is computed within 
the framework of the Foldy-Lax equation, refer to 13j, |24|, |27|. Then, a white Gaussian 



noise with 20dB signal-to-noise ratio (SNR) is added to the unperturbed data in order to 
investigate the robustness of the proposed algorithm while at least partly alleviating inverse 
crime. In order to obtain the number of nonzero singular values Mf at each frequency Uf, 
a 0.01-threshold scheme (choosing first j singular values Sj(uf) such that ~4^ry > 0.01) is 
adopted. A more detailed discussion of thresholding can be found in 23|, |24j, |25| . As for the 
step size of the search points x 6 tt, it is taken of the order of 0.02. 

A is observed from the numerical experiments led in [6|, similar results were obtained 
for both cases e + > and e + = e_ (permeability, and both contrast cases also). Thus, two 
different situations of interest are henceforth considered. 
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4-1. Permittivity contrast case: St ^ £- and \it — A 4 - = /-*+ 
4.1.1. e + > £- case 

In this example, we choose the values e + = 5 and e_ = 4. Thus, the number, iV + , of 
transmitted waves is smaller than the one, N, of directions of incidence. Based on the test 
configuration in Table [TJ imaging results are depicted in Fig. [3j A good imaging resolution 
appears to be reached when the thin inclusion has a small constant curvature, e.g., Fig. 



3(a) , but this resolution is poorer when it is of large curvature, e.g., Fig. 3(b) Nevertheless, 
we can still conclude that a thin inclusion has been successfully retrieved. 

The proposed algorithm could be applied directly to a target consisting of several well- 
separated thin inclusions with same thickness h. Let us skip the derivation herein and 
only provide imaging results. Those are depicted in Fig. [3] with Ti and T 2 , the respective 
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inclusion 


N 


N + 


F 


frequency 


range 


search domain 




36 


28 


32 




0.2 J 


ft 


= [-1,1] x [-3, 


-1] 


r 2 


42 


36 


36 




2tt 
0.2 J 


ft 


= [-1,1] x [-3, 


-1] 


r 

1 M 


50 


40 


42 


^[fe 


2tt 
0.1 J 


ft 


= [-1,1] x [-3, 


-1] 



Table 2: (Permeability contrast case) Test configuration for Ti, T2 and T M . 



permittivities being e\ and 82- We notice that, if an inclusion has a much smaller value 
of permittivity than the other, it does not significantly affect the scattering matrix (as 
expected) and consequently, it cannot be retrieved via the proposed algorithm, refer to Fig. 



3(d) (see [24], Section 4.5] also). 



4-1-2. e + < e_ case 

In contrast with the previous example, let us consider a more practical situation where 
the upper half-space (air) is the least refractive. In this case, we let e + — 1 and e_ = 3. 
Now, the number of directions of transmitted waves, iV + , is the same as the one of directions 
of incidence, N. 

Figure H] shows the result based on the test configuration in Table [TJ Similarly with the 
previous example, a thin inclusion of small constant curvature is well retrieved. However, 
poor results are observed when the inclusion is of large curvature. It is interesting to observe 
that the location of the end-points of Ti and r 2 is well identified however. That is, connecting 
them by a straight line, it should provide a good initial guess for an iterative solution 



algorithm. Similarly with the Fig. 3(d) whenever an inclusion has a significantly smaller 



value of permittivity than the other, it cannot be successfully retrieved via the proposed 



algorithm, refer to Fig. 4(d) 



4-2. Permeability contrast case: Et = £- = £+ and fir 7^ A*- 

Similarly with the previous case, two different situations of interest are considered based 
on the test configuration in Table [2] 

4-2.1. fi + > /u_ case 

In this example, we choose the values //+ = 5 and /x_ = 4. Thus, we have a number, 
iV + , of transmitted waves which is smaller than the one, N, of directions of incidence, b 
is set to cb = (0, 1) (see Remark l3.ip . Based on the test configuration in Table [21 imaging 
results are depicted in Fig. [5J Similarly with the permittivity contrast case, a good imaging 
resolution is attained when the thin inclusion has a small constant curvature but a poorer 
resolution when it is of large curvature. In the result corresponding to multiple inclusions, 
similarly with what happens in the permittivity contrast case, we easily notice that if an 
inclusion has a much smaller value of permeability than the other, it cannot be retrieved via 



the proposed algorithm, refer to Fig. 5(d) 
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11 




12 






4-2.2. fx + < [X- case 

In contrast with the previous example, let us consider that the upper half space is the 
least refractive. In this case, we set /i + = 1 and //_ = 3. 

Results are shown in Fig. E] based on the test configuration in Table [2J Similarly with 
the e + < E- case case (Fig. @J, the location of the end-points of the thin inclusion is well 
identified. As for the imaging of multiple inclusions, tough we cannot reach a proper result, 
a good initial guess is made possible. 

4-3. Permittivity and permeability contrast case: Et ^ £- and fir ^ Ai_ 

Similarly with the previous cases, two different situations of interest are considered based 
on the test configuration given in Table [3j 
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inclusion 


N 


N + 


F 


frequency 


range 


search domain 




40 


32 


36 




0.2 J 


ft 


= [-1,1] x [-3, 


-1] 


r 2 


48 


40 


40 




2tt 
0.2 J 


ft 


= [-1,1] x [-3, 


-1] 


r 

1 M 


60 


50 


48 


^[fe 


2tt 
0.1 J 


ft 


= [-1,1] x [-3, 


-1] 



Table 3: (Both permittivity and permeability contrast case) Test configuration for Ti, T2 and T M . 



4-3.1. e + > 8- and /i+ > //_ case 

In this case, the values e + = 5, /i + = 5, e_ = 4 and //_ = 4 are chosen, i.e., c is set to 
c = (1, b) = (1,0, 1) (see Remark [3.1ft . Based on the test configuration in Table [3] imaging 
results are depicted in Fig. [0 In opposition with the previous cases, a poorer resolution is 
now achieved. However, at least, we can approximate the shape of the thin inclusions unless 
one has a much smaller value of permittivity and permeability than the other. 

4-3.2. e+ < £_ and /i + < //_ case 

Now, let us consider that the upper half space is the least refractive. In this case, we let 
£ + = 1, fx + — 1, £_ = 3 and //_ = 3. Based on the test configuration in Table EJ imaging 
results are depicted in Figs. [8] Although poor results are found when the inclusion is of 
arbitrary shape, opposite to the previous examples, we still recognize a part of the thin 
inclusion r 2 (neighborhood of point (0.3, —2.5)) has a much smaller value of permittivity 



and permeability than the other, refer to Fig. 8(d) 



5. Conclusion 

A non-iterative imaging algorithm operated at several time-harmonic frequencies has 
been proposed in order to determine the unknown support of a thin penetrable electromag- 
netic inclusion buried within a half space. The main idea comes from the fact that the 
collected Multi-Static Response (MSR) matrix data can be modeled via a rigorous asymp- 
totic formulation for a limited range of incident and observation directions. 

Thanks to various numerical simulations, it is shown that the proposed algorithm is 
both effective and robust with respect to random noise. Moreover, it can be easily applied 
to multiple inclusions. Nevertheless, some improvements are still required, e.g., when the 
sought inclusion is of large curvature or has a significantly smaller value of permittivity (or 
permeability) than, say, another nearby. In addition, in order to achieve a proper imaging 
of arbitrarily shaped thin inclusions in the case of permeable, or dielectric and permeable 
material, the choice of the testing vector functions b and c still requires further mathematical 
investigation. 

It is worthwhile emphasizing that such results are obtained at low computational cost 
but that they do not guarantee complete shaping of the inclusions. However, they could 
provide a good initial guess of a level-set evolution [l], QjJ 26] or of any other standard 



iterative algorithm [14 . 
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-1 -0.5 0.5 1 -1 -0.5 0.5 1 

(c) r M with £i = Ei = fix = /12 = 5 (d) r M with e\ = fii = 10 and £2 = A*2 = 5 

Figure 7: (e + > £_ and n + > /i_ case) Map of W(x) for Ti, T 2 and r M . 
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(a) Ti with £i = /ii = 5 
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0.2 
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0.6 
0.5 
0.4 
0.3 
0.2 
0.1 



-0.5 0.5 1 

(c) T M with £i = £2 = Mi = /"2 = 5 



-1.5 



-2.5 



-0.5 0.5 1 

(b) T 2 with e 2 = [i2 = 5 



n 



J 



0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




-0.5 0.5 

(d) r M with £i = fii = 10 and £2 = A*2 = 5 



Figure 8: (e + < £_ and fi + < /i_ case) Map of W{x) for T\, T 2 and Tn 
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Finally, we have been considering a two-dimensional problem. The strategy which is 
suggested, e.g., mathematical treatment of the asymptotic formula, imaging algorithm, etc., 
could be extended to the three-dimensional problem, refer to 0, Est for related work. 
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